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Abstract 

We consider here variational solutions in the Hartree-Fock approximation upon 
breaking time reversal and axial symmetries. When decomposed on axial harmonic 
oscillator functions, the corresponding single particle triaxial eigenstates as func- 
tions of the usual cylindrical coordinates (r, 9, z) are evaluated on a mesh in r and 
z to be integrated within Gauss-Hermite and Gauss-Laguerre approaches and as 
Fourier decompositions in the angular variable 9. Using an effective interaction of 
the Skyrme type, the Hartree-Fock hamiltonian is also obtained as a Fourier se- 
ries allowing a two dimensional calculation of its matrix elements. This particular 
choice is shown to lead in most cases to shorter computation times compared to the 
usual decomposition on triaxial harmonic oscillator states. We apply this method 
to the case of the semi-quantal approach of large amplitude collective motion cor- 
responding to a generalized routhian formalism and present results in the A = 150 
superdeformed region for the coupling of global rotation and intrinsic vortical modes 
in what is known after Chandrasekhar as the ^-ellipsoid coupling case. 



1 Introduction 

In microscopic calculations within the Hartree-Fock approximation, one-body 
eigenstates wave functions can be either evaluated numerically on a mesh or 
decomposed on some truncated set of orthogonal functions. In the former case, 
it is necessary to solve the static Schrodinger-like integro-differential or in the 
Skyrme case partial-derivative equations which can prove to be a long process, 
whereas in the latter, one has to compute matrix elements of the Hartree- 
Fock hamiltonian on the set and then diagonalize the obtained matrix. This 
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second method is known to be rather efficient in terms of computation time. 
However, as one is bound to use only a truncated basis, the obtained results are 
dependent on the basis parameters and size. Nevertheless the energy difference 
between the truncated solution and the exact one (which would correspond to 
an infinite basis and can be approached by the former method) may in some 
cases be evaluated (see e.g. [1]). 

In order to describe a compact system of particles interacting through mostly 
attractive two body forces, one often chooses as a basis set, eigenfunctions 
of a deformed harmonic oscillator single particle hamiltonian. One truncates 
it for instance [1] by only retaining the eigenstates whose eigenenergies are 
lower than some truncation energy reference. The obtained solutions thus 
depend on both the truncation energy and the deformation parameters of the 
harmonic oscillator basis. For a given truncation, the deformation parameters 
to be used are obtained by a minimization of the energy with respect to these 
parameters. When describing axially symmetric nuclei, it is appropriate to 
use axially symmetric harmonic oscillator wave functions, all the calculations 
being then performed in two dimensions. To the best of our knowledge, when 
projecting triaxially shaped solutions on a single particle basis, only triaxial 
harmonic oscillator wave functions have been used [2-4]. We present here 
an alternative way to describe triaxial shapes by using an axially symmetrical 
basis. Axially symmetric wave functions are eigenstates of the projection of the 
angular momentum on the quantization axis. By mixing states with different 
eigenvalues for this operator, one can obtain triaxial wave functions. As a 
matter of fact, in the position representation, those functions are expressed 
for any point of the (r, z) plane as a complex Fourier series in the angular 
variable 9. 

Using, as we will do here, an effective interaction of the Skyrme type (namely 
in the SkM* parametrisation [5]), the total energy of the nucleus is obtained 
by integrating an energy density functional which, as well as the potentials 
of the Hartree-Fock hamiltonian, is expressed in terms of sums, products and 
integer and non-integer powers of some local densities. With our choice of 
axially symmetric basis states as we will show, those densities, and hence the 
Hartree-Fock potentials and energy functional, are obtained as real Fourier 
series in the angular variable. By integrating the energy functional over the 
whole space, only the first term of the Fourier series (which remains constant as 
a function of 9) is non vanishing. This component is calculated for each point in 
the (r, z) plane and then integrated in two dimensions. Using selection rules, 
the calculation of the Hartree-Fock hamiltonian matrix elements can also 
be performed through a single two-dimensional integral whereas when using 
triaxial wave functions, all these integrations would be three-dimensional. Our 
alternative method thus is expected to lead to shorter computation times 
provided that the number of axially symmetric basis wavefunctions necessary 
for an accurate description of the triaxial solution is not too large. 
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Even for triaxial shapes, high order Fourier components are exactly zero and 
the number of non vanishing components is determined by the truncation 
parameter. In practical cases, one can impose a cutoff in the maximal order of 
the retained components, due to the observed fast convergence of the numerical 
results with respect to this maximal order. For axially symmetric nuclei, only 
the first component of the densities are non vanishing. We are thus able to 
describe both axially symmetric and triaxial shapes within the same code by 
simply adjusting the maximal order of the Fourier decompositions. For largely 
triaxial shapes we can expect relatively important truncations effects within 
our method, which will lead to a reduced competitivity of our approach with 
respect to those using triaxial basis. 

This calculation method has been applied to physical situations of triaxial 
shapes combined with time-reversal symmetry breaking as obtained within 
the well known cranking (or routhian) approximation or within approaches 
corresponding to a generalization of the latter. This generalization which will 
be called below the generalized routhian formalism has been briefly sketched in 
a previous paper [6] and will be presented in more details in Ref. [7]. It allows 
us to describe the dynamics of a class of collective modes defined by a velocity 
field in a semi-quantal approach a la WKB. It reduces to the addition of a 
time-odd constraint — (3 ■ p to the static hamiltonian, where (3 is the classical 
collective velocity field under study and p is the single particle momentum 
operator. The possible choices for the /3-fields have been restricted by imposing 
the routhian eigenstates to be also eigenstates of the parity and of the signature 
operator with respect to the first axis. 

This paper will be organized as follows. In Section 2 we will present our formal- 
ism with some calculation and computation details. Section 3 will be devoted 
to discussions about truncation effects and the computation time we expect 
to gain and will end with the comparison of our results with the literature. We 
will then turn to the generalized routhian formalism in the .S-type ellipsoid 
case and present some preliminary results in section 4 before concluding and 
presenting perspectives of future work in section 5. 



2 Hartree— Fock— Skyrme generalized routhian 

2.1 One and N-body routhians 

Within classical mechanics, the canonical form of the equations of motion 
is conserved under a canonical transformation of the dynamical variables by 
adding to the hamiltonian the time derivative of a generating function. In the 
case of a local transformation (i.e. involving only the coordinates), this time 
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derivative writes as 



where (3 is the collective velocity related to the transformation. In the quantum 
mechanical case, it can be shown [7] using a certain class of unitary transforma- 
tions that the density matrix solution of the Hartree-Fock equation involving 
the generalized routhian 

R = H-^((3-p + p-(3) = H-(3-p+%iv(3. (2) 

in lieu of the hamiltonian H exhibits in the laboratory frame a dynamical con- 
tent which corresponds to the collective velocity field (3. When spin degrees of 
freedom are present, it is necessary to also consider the spin-rotation collective 
mode. Considering a vortical velocity field, the routhian then writes 

R = H-(3-p-^n-cr, (3) 

where a is the vector whose components are the three Pauli matrices. As an 
example, in the case of global rotation where the collective field is 

ft* = n x r, (4) 

the routhian is expressed by 

R = H-(3 TOt -p-^n-* = H-n-(£ + s), (5) 
which is the well known cranking approximation. 

Using Skyrme interaction, the routhian expectation value for a Slater deter- 
minant can be written as the integral over the whole space 



(R}= J (x(r)-hf3-j-^n-p^j dr 



(6) 



where % is the energy density functional, j is the current density and p is the 
spin- vector density whose expressions are given in Appendix A.l. 

Minimizing this expectation value with respect to all single particle wave func- 
tions, one obtains the following expression for the one-body routhian (see 
Ref. [8] for instance) 



h q = U g - ^- Vf g ■ V + j (a q ■ V + V • a q ) 



-n (s q + iwy;° x v) • <r, (7) 
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where q stands for the considered charge state. The expressions of the various 
form factors entering this routhian are given in Appendix A. 2. 

2.2 Symmetries 

As said in the introduction, we have chosen here to use axially symmetric 
harmonic oscillator wave functions for our basis states. Such wave functions 
have been widely used [9] to describe axially symmetric variational solutions. 
For the single particle hamiltonian, they entail symmetries with respect to 
the parity and the third component of the angular momentum, and can be 
written in the coordinates representation in terms of normalized Hermite and 
associated Laguerre polynomials H nz and (see Ref. [10]) as 

1/2 e iA ^ |A|/2 ^(04 A J(^), (8) 

where /i represents the set (n z n r A) of quantum numbers, f3 z and j3± are the 
usual oscillator constants [9] which are generally related to the parameters 
A) — (PzPi) 1 ^ 3 and q = (f3±/P z ) 2 and £ and r\ are given in terms of the 
cylindrical coordinates r and z by 

£ = z(3 z , V = r 2 (3 2 ± . (9) 

For notational convenience, we will omit in the following the £ and r\ depen- 
dence of the polynomials. 

Whenever the time-reversal symmetry is broken, the third component of the 
angular momentum can no longer be chosen as a symmetry. It is well known 
however, that the Hartree-Fock hamiltonian obtained from the Skyrme in- 
teraction, the parity and for instance the first component of the signature 
operator [11] defined as 

S 1 = U 2 U 3 a u (10) 

(where II; is the reflection operator in the % direction, and o\ is the usual Pauli 
matrix) form a set of commuting operators. Imposing that the routhian also 
commutes with the two last symmetry operators limits the possible choice for 
f3. As an example, if we restrict ourselves to first order polynomials in the 
coordinates, the vortical velocity field given by its components 

A = o, 

= ax 3 , (11) 
As = bx 2 



yv( r ) 



7Y 
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is the most general one that fulfills these commutation requirements. This an- 
alytical form corresponds to the so-called Riemann 5"-type ellipsoid solution 
presented by Chandrasekhar [12] in the context of fluid dynamics, whose ap- 
plication to nuclear physics has been discussed in particular in [6,13]. This field 
can then be rewritten as the sum of a global rotation term of angular velocity 
Q and an intrinsic vorticity term of angular velocity uj both perpendicular to 
the first axis in the form 

A = o, 

fo = -(n + u/q)x 3 , (12) 
/3 3 = (Q + uq)x 2 , 

where q is the as/a 2 axis ratio of the ellipsoid approximating the nuclear shape. 
The coupling of the rotation perpendicular to the first axis with the spins is 
taken into account by a term proportional to fieri in the routhian which still 
commutes with the signature operator. 

Since the routhian single-particle eigenstates must have the same symmetries 
as the routhian, it is appropriate that the basis states are also eigenstates of 
the parity and signature operator. The action of the parity P and signature 
Si operator on axial harmonic oscillator states is given in usual notations by 



P\n z n r kT) = (-) nz+A |n z n r AS), 
Si \n z n r AS) = (— ) nz \n z n r — A — E) , 



(13) 
(14) 



from where we can see that Si 2 is the identity operator, and thus the eigen- 
values of Si are ±1. We then obtain eigenstates of the parity and signature in 
the form 



\^) = ^\H) + s(-rm-±) 



s = ±1, 



(15) 



where \i {p) represent the set (n z , n r , A) [(n z , n r , —A)] of quantum numbers, 
and we have 



P\H8) = (-r +A |/i*>, Si\pis) = s\pis). 



(16) 



The single-particle eigenstates obtained after minimizing the routhian are 
written as 

l*> = ( 17 ) 

where the sum runs over basis states having the same eigenvalue for the two 
symmetry operators. In the coordinates representation, these eigenstates can 
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be decomposed in two components with spin | and — | as 

^(r) = *if(r)||>+^(r)|-i>, 
and one deduces from equations (8) and (17) 



(18) 



*£(r) = 



PS 

2tt 



nl/2 



E(±) n ^; e±iA ^ |A|/2 ^4 A j- (19) 



We are now able to calculate the expressions of the various densities of Ap- 
pendix A.l. The density p of equation(A.7) for instance rewrites with the 
notations of (18) as 



(20) 



since s 2 — 1. Now inserting the expression (19) and rearranging the terms we 
get 



P = 



2^" \ k 



e -(^) E fc^C* ) ^(|AHA'l)/2^ L jA| ffB;L IA'l 



jgi(A'-A)0 + (_)n,+T»; e i(A-A')^ _ (21) 

The sum over /c only involves the components C k of (17) and is nothing but 
the matrix element of the density matrix, that is 

ECf^/W- (22) 

k 

In the general case this matrix is hermitian, but it can be shown [14] that the 
choice made for the basis states makes it a real symmetrical matrix. Using 
this property, the imaginary parts of the pp' and p'p terms in equation (21) 
cancel each other, and we can write 



P = 



P Z P 2 



p -(f+l)Vfl ,„(|A|+|A'|)/2it t\M H ,t) A 



|A'| 



2tt 



cos[(A - A')6] (l + (-)"*+<) . (23) 

Finally, since the basis states \p) and \p') have the same eigenvalue for the 
parity operator we can deduce from (16) that 

(_)»,+n; = (-) A ~ A ', (24) 

and the pp' term of the sum (23) vanishes if the difference A — A' is an odd 
integer. We then write at last the density p as 
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Table 1 

Definition of the non-vanishing components of the Fourier series representing the 
various local densities in use; their parity tt z with respect to z is also given. 



Densities 


function 


7T 2 


p, r, V 2 p, V • J 


cos 2p9 


+ 


j r , cur hp 


sin[(2p + 1)0] 




j 6 , curl e p, curl 2 j, p z 


cos[(2p + 1)9] 




j z , curl 2 p, curl^j, p e 


sin[(2p + 1)6»] 


+ 


curl r j, p r 


cos[(2p + 1)9] 


+ 



?2 

c2_ 



P = ^e"^ £ P^ m+m)/2 ^-v cos[(A - A')0] 
n,n' 

H nz H n ,Ll A }L[f. (25) 
with the definition 



^ 2p i 1 if A — A' is an even integer, 
if A — A' is an odd integer. 



The calculations of the other local densities entering the one-body routhian 
are presented in Appendix B.l. It is clear from equation (25) and from the 
results of this Appendix that all the densities are obtained as real Fourier 
series in the angular variable 9. As some of the Fourier components of these 
densities are vanishing, we recall in Table 1 the non-vanishing components 
together with their parity with respect to the z variable. The maximal order 
of the non- vanishing components is 2A max + 1 where A max is the higher A value 
among the basis states which is equal to in the case of a spherical basis 
and increases for deformed basis. 

The form factors of the one-body routhian are expressed in Appendix A. 2 as 
sums and products of these densities, non-integer powers of p and the Coulomb 
potential. We show in Appendix B.2 how the Coulomb potential can be ob- 
tained as a Fourier series in the 9 variable, and the Fourier decompositions 
of the non-integer powers of p are obtained through weighted integrals over 9 
computed with a 48 points Gauss-Legendre quadrature formula. The Fourier 
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Table 2 

Same as Table 1, but for the one-body routhian form factors. 



Form factors 


function 


z 


U, f, V so 


cos 2p6 


+ 


a r 


sin(2p + 1)6* 




a e , S z 


cos(2p + 1)0 




a z , S e 


sin(2p + 1)6* 


+ 




cos(2p + 1)6 


+ 



decomposition of (3 is readily obtained from equation (12) as 

P r = zip + ujq) sin#, 

Pe = -z(n + ujq)cos6, (27) 
P z = —riVL + uj/q)sm6. 

Using the well known relations expressing products of trigonometrical func- 
tions as sums of trigonometrical functions of the sum and difference of their 
arguments, it is then straightforward to obtain the routhian form factors as 
Fourier series, whose no n- vanishing components are given in Table 2. 

Having expressed the routhian form factors as Fourier series in the angular 
variable allows us to calculate the routhian matrix elements as two-dimensional 
integrals, using simple selection rules. Indeed, if we write for instance, referring 
to Table 2, the scalar potential U as 

u(v,e,o = £ 008(2^)^(17,0, (28) 

n>0 

the part of the matrix element involving this form factor is spin diagonal and 
writes 



(fis\U\n's) = \ £>| cos(2n#) U (2n) \fi') + i-) nz+n '*(fl\ cos(2n0) C/ (2n V), 



n>0 



(29) 



readily performing the spin part of the scalar product. For a given n, the 
angular part of the first term on the right-hand side is simply proportional to 



r2n 

Jo 



; , , . , ,n 1 if /// / A A' . . , 

V(1 + <^ A _ A ,|) ifm=|A-A'|, V ; 
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and is identical to the angular part of the second term up to a complex con- 
jugation. It can then be easily seen that equation (29) reduces to 

(HtVs) = H\ cos (l A - A '^) f/(|A ~ A ' 1 Vi>> ( 31 ) 

which writes in analytical form 

P OO P oo 

(fis\U\fi's) = (1 + 5°_ A ,) y dry y d^^V^^'l)/ 2 

#" A %J Bl «, (32) 

where we have restricted the £ integration to positive values, since the in- 
tegrand is an even function of £. The calculation of the other terms of the 
routhian matrix elements also reduces to two-dimensional integrals and is 
presented in Appendix C. These integrations are performed through 10 points 
Gauss-Hermite and Gauss-Laguerre integrations (for the Gauss-Hermite in- 
tegration, we use in fact a 20 points method which reduces to 10 points due 
to n 3 symmetry properties). We therefore only need to compute the values 
of these densities and routhian form factors on the mesh points of a 10 by 10 
grid. 

2.3 Calculation of some observables 

Within the Skyrme-Hartree-Fock formalism, the calculation of observables is 
generally reduced to a three-dimensional integral as is obviously the case for 
the routhian expectation value in equation (6). Expressing this integral as the 
integral of a Fourier series, only the first component needs to be integrated - 
which is done through a two-dimensional integral — since the integrals of the 
others are vanishing. For instance, the root mean square radius of the mass 
distribution which is given by 

RLs = \\ r 2 P (r) dr (33) 

where A is the total number of nucleons, reduces to the two-dimensional in- 
tegral 

Aqr roo roo 

RLs = ^J o rdrj o dz(r 2 + z 2 )p^(r,z), (34) 

since r 2 does not depend upon 6. The quadrupole operator expectation value 
is easily obtained as 

poo roo 

(Q ) = 4tt / rdr dz(2z 2 - r 2 )p<®(r,z), (35) 
Jo Jo 
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similarly to what has been done for the root mean square radius calculation. 
The a 2 semi-axis length is given up to some normalization factors as the square 
root of (x 2 2 ) written as 

(x 2 2 ) = jjr 2 sin 2 6p{r) dr. (36) 

The integrand easily writes as a Fourier series upon developing sin 2 6, and 
using the results of Table 1. Retaining only its first component, we get 

(a; 2 2 ) = ^lj\drj^ dz r 2 (p(°> (r, z) - \p {2) (r, z)) . (37) 

And the expectation value of the non-axial quadrupole moment Q 22 is obvi- 
ously 

(Q 22 ) = A(x 2 2 - x ± 2 ) = -2n rdr dz r 2 ' p {2 \r , z) . (38) 

JO JO 



The first component of the total angular momentum is obtained from the 
current and spin-vector densities as 

{L 1 /h) = ((rxj + p)-e 1 ), (39) 

where e x is the unit vector in the first direction. It is expanded as 

(Li/K) — J rj z sin 9 — z{je cos 6 + j r sin 6) + \{p r cos 6 — p$ sin 6) dr, (40) 



and finally, referring to Table 1 we get, 



9-7T roo roo 

(LJh) = — rdr dz 
A jo jo 



W+j^ + WP-pP)]- (41) 



Taking into account the total angular momentum quantization rule for an even 
number of nucleons, we will only retain solutions that satisfy 



{L./h) 2 = 1(1+1), 



(42) 



where / must be an even integer. The dynamical moment of inertia is cal- 
culated numerically as the first derivative of this angular momentum with 
respect to the angular velocity, namely 



CJ(2) _ ?L 



(43) 



which, as mentioned by Gall et al. [15], provides a better numerical stability 
than the second derivative of the energy with respect to Q. 
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3 Numerical tests 



3.1 Convergence and computation time 

Describing as we do here the single-particle states as a decomposition on a 
truncated harmonic oscillator basis, it is then crucial to study carefully the 
convergence of the solution with respect to the truncation. We are using the 
usual deformed truncation scheme, namely considering only as basis states 
those having an energy lower than what is obtained for a state with N quanta 
of the equivalent spherical oscillator [9]. Discussions on this problem in the 
static case can be found in the literature [1] . We will thus concentrate here on 
the rotational case. 

We have studied as an example rotating solutions of the nuclei 80 Sr. In Table 3 
we present some properties of the I = 20 h solution for No varying between 6 
and 14 major shells. The deformation parameters were obtained from an op- 
timization of the static solution. One sees that the routhian expectation value 
which is converged up to 10 eV (for a given basis size) vary by less than 0.23 
percent as one adds two more shells. We have checked that the convergence of 
the various terms entering the routhian (eqs. A.2-A.6) also converge with such 
a good accuracy, with the noticeable exception of the time-odd contribution 
which is roughly 700 times smaller than the total routhian and varies by 5 
percent between 8 and 10 shells. Now looking at the deformation properties, 
we see a rather nice convergence of the quadrupole moment above Nq = 8. 
The variations of the non-axial quadrupole moment are of the same order in 
absolute values. 

In the lower part of Table 3 we can see that the angular velocity varies by less 
Table 3 

Values of some observables in the I = 20 h solution of the nuclei 80 Sr as a function of 
Nq (see Subsection 2.3 for calculation details). We used the parameters (3q = 0.534 
and q = 1.2658. 





6 


8 


10 


12 


14 


(R) (MeV) 


-686.69 


-688.07 


-688.75 


-690.28 


-690.61 


(Qo) (b) 


6.04 


6.10 


6.12 


6.13 


6.10 


(Q22) (b) 


0.32 


0.44 


0.47 


0.41 


0.43 


Q (MeV) 


0.831 


0.831 


0.831 


0.836 


0.832 


3?( 2 ) (^.MeV" 1 ) 


24.46 


24.06 


23.96 


23.93 


23.91 
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Table 4 

Computation time in seconds (measured on a HP9000/780 workstation), of one 
Hartree-Fock iteration within various formalisms as a function of Nq. In the triaxial 
case, the first iteration is referred by W and the following by ( 2 ) . 





6 


8 


10 


12 


14 


Vautherin's formalism 


0.24 


0.59 


1.46 


3.07 


5.98 


"axial" case(parity,signature) 


0.57 


1.87 


7.31 


32.8 


108 


triaxial case ^ 


9.65 


15.6 


28.8 


68.5 


169 


triaxial case ( 2 ) 


1.49 


4.69 


15.2 


52.1 


150 



than 0.5 percent as two shells are added. The convergence of the moment 
of inertia is even better, since its variation rate is constantly decreasing. This 
good convergence of the moment of inertia can be easily interpreted. Indeed, 
the moment of inertia is a differential quantity, and it is reasonable to expect 
small truncation effects for it in the near vicinity of a given solution. 

Let us now discuss the computation time required by our formalism. First of 
all, we present in Table 4 the time needed to perform an iteration in the axial 
and the triaxial case, that is in the limiting case where only the first Fourier 
components are retained and in the case where all the non-vanishing compo- 
nents are used. In the triaxial case, we give the values for the first iteration 
where some initialization parameters are computed (see Appendix B.2) and 
for the following iterations which are shorter. We also present for comparison 
the computation time used by the static (and axial) formalism of Vautherin 
[9]. In static cases, the ratio between our formalism and Vautherin's is some- 
how small for low N values, but almost reaches 20 for iVo = 14. This is due to 
our choice of the parity and signature symmetries for the single-particle states 
which only split each parity block into two, in contrast to Vautherin's case. 
The time reversal symmetry breaking also doubles the size of the basis. On 
the contrary, excluding the initializations performed during the first iteration, 
the ratio between the triaxial and axial case in our formalism always decreases 
and stays relatively small. 

An almost direct comparison of our code could be performed with the one built 
by Dobaczewski and Dudek [16] using a triaxial harmonic oscillator basis for 
the 152 Dy nucleus with the SkM* effective force. The basis size is defined by M 
of the order of 300 in their notation. This corresponds roughly to a deformation 
dependent truncation with N = 10 in our notation. These authors quote a 
time of 9 seconds on a CRAY C90 for the vectorised version. They provide a 
factor of two to three to translate into the time used on e.g. an IBM RS/6000 
station for a non vectorised version. The corresponding time in our calculation 
of 150 Gd with roughly the same basis size and after the first iteration, would 



13 



Table 5 

Computation time in seconds of various parts of our formalism during the first iter- 
ation (measured on a HP9000/780 workstation). The last line gives an expectation 
of the time required by a fully triaxial formalism using triaxial wave functions. 





6 


8 


10 


12 


14 


Diagonalization 


0.08 


0.48 


3.19 


22.4 


83.2 


Coulomb potential 


8.17 


10.9 


13.6 


16.4 


19.1 


Matrix elements 












and densities 


1.05 


3.63 


10.6 


25.8 


56.4 


Fully triaxial 


13.7 


40.1 


113 


289 


661 



be of less than 20 seconds which is of the order or slightly better than the 
quoted time. 

Another way to estimate the advantages, or at least the competitiveness of 
our approach consists in extrapolating from our code the computation time re- 
quired by triaxial basis calculations. For instance, in our case the routhian ma- 
trix element calculation involves two-dimensional integrations, where it would 
involve three-dimensional integrations in the fully triaxial case. Using 10 points 
Gauss-like integration method we can expect the matrix elements to be com- 
puted 10 times faster within our formalism. The same factor is expected for the 
densities calculation. In our formalism, each of the N non-vanishing Fourier 
components of the Coulomb potential are calculated on a 10 by 10 mesh in 
the (r, z) plane performing a 10 by 10 integration in the (r', z') plane and a 48 
points integration in the angular variable 9. In the fully triaxial case, it has to 
be computed for each of the points in the (x, y, z) mesh, each time performing 
a three-dimensional integral with 10 points. We thus expect our formalism to 
be roughly a factor Nq/2 slower than the fully triaxial case during the first 
iteration. The routhian diagonalization is also time-consuming but it should 
take the same amount of time within the two formalisms. Finally it seems that 
the computation of the energy and the routhian form factors should be slower 
in our formalism, but the time required for those calculations is negligible as 
compared to the overall execution time of an iteration. 

All these remarks allow us to calculate an "expected" computation time for 
the fully triaxial formalism. We present this expectation in Table 5 together 
with the computation time of the discussed parts of our formalism. Comparing 
the last line of Tables 4 and 5, we see that our formalism is expected to be 4 
times faster for Nq = 14, and this factor is likely to be even larger for lower 
N values, reaching a factor of 9 for N = 6. These factors result of course 
only from a guess, and should be confirmed by existing fully triaxial codes. 
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Of practical interest also, would be a comparison between the computation 
times needed in our case and in other competing and totally different ap- 
proaches to describe heavy nuclei within the same range of deformations. 
For the same 150 Gd nucleus at the same deformation (ground-state (3 value, 
7 = 15°) and for the same basis size (N = 8) and iteration number (30), the 
Bruyeres le Chatel code (discarding the time needed to treat pairing correla- 
tions) would take [17] about 30 percent more time than our code. Of course, 
the Gogny effective force which is used there is completely different from our 
Skyrme force. Another point of comparison could be the 3D calculations of 
solutions having the same symmetries for slightly heavier deformed nuclei 
(A ~ 190) using the imaginary-time step technique where the computing time 
for a similar degree of convergence would be slightly larger [18] than ours. 
However it should be stressed that in their case the time used to treat pairing 
correlations has not been discarded. 

It then appears that the orders of magnitude of the computing times for all 
considered codes, are identical with a slight advantage of our method in com- 
parable approaches. It could be noted finally that the rather recent character of 
the present version of our code should reasonably leave some place for further 
optimization. 



3.2 Results of superdeformed 150 Gd 

We have chosen to test the physical results obtained within our formalism by 
comparison with the Hartree-Fock results of Bonche et al. [19] for the yrast 
band of the nucleus 150 Gd as calculated with the SkM* parameterization. 
Their formalism is based upon the computation of the one-body eigenstates 
on a three-dimensional mesh by solving the Schrodinger-like Hartree-Fock 
equation and thus are free from truncation effects but of course are contingent 
upon mesh size and related approximate numerical derivatives. The calcula- 
tions presented here have been performed using 10 major shells. Though we 
should take into account at least the 21 lower order Fourier components of 
the densities, we have checked that the results are unchanged if we only use 7 
Fourier components. This is due to the very little triaxiality of the solutions 
whose maximal triaxial axis ratio is only 1.02 (corresponding to 7 < 1°) and 
is reached at the end of the band. 

The basis parameters optimization has been performed in several steps. We 
have obtained a first set of parameters from an optimization of a static solution 
constrained to the experimental value of the charge quadrupole moment of 
the yrast band. We have then minimized the routhian expectation values for 
rotating solutions with angular momentum ranging from 10 to 80 h with steps 
of 5 h. Interpolating between these points, we obtained optimized parameters 
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Fig. 1. Dynamical moment of inertia of the Gd yrast band as a function of the 
angular velocity. The two curves at the bottom of the figure are the contribution of 
neutrons (upper curve) and protons (lower curve) to the moment of inertia. 

with steps of 2 ft. 

We finally minimized the routhian for each of these quantized values of the 
angular momentum reaching an overall precision on the minimum of less than 
0.1 keV. All these optimizations have been performed through successively 
minimizing quadratic fits of the routhian in the /3q and q variables. We only 
performed calculations for angular momentum above 10 ft because we did 
not want here to describe the transition from superdeformation to normal 
deformation. 

We show in Figs. 1 and 2 respectively the Dynamical moment of inertia and 
the charge quadrupole moment (calculated by replacing the matter density 
p in eq. (35) with the proton matter density p p ). While our values of the 
dynamical moment of inertia is in perfect agreement with the one depicted 
in [19], we obtain a charge quadrupole moment which is slightly larger by 0.2 
e 2 b. This effect cannot be related to truncation effects since we have checked 
that the same quadrupole moments are obtained with N = 14. However, the 
variation of this moment over the yrast band is quantitatively well reproduced 
by our calculations. 

We finally plot in Fig. 3 the neutron and proton single-particle routhians. We 
are globally in good agreement with the results of Ref. [19] (see for instance 
the proton gap for Z = 66 between the [651] and [411] states and the intruder 
orbital in the proton spectra appearing on top of the spectra above hQ = 
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Fig. 2. Charge quadrupole moment of the 150 Gd yrast band as a function of the 
angular velocity. 

0.5 MeV) but we observe some small differences. Among the less important are 
the labels of the eigenstates which are directly obtained in our formalism from 
the decompositions on the axial basis. In the neutron spectra, our [521] 3/2 
eigenstate (with a squared overlap of 60 percent) is labeled [512]3/2 by Bonche 
et al., and in the proton spectra, their [532] 5/2 eigenstate is found to be in our 
formalism [303]5/2 (with a squared overlap of 81 percent) while the [523]5/2 
state lies 500 keV below. 

Comparing more closely the two results, the level spacings are seen to slightly 
differ, leading to different values for some level crossings. In some cases, such 
small effects which we believe to be related to numerical methods in both cases 
can appear to be relevant, changing the properties of a given band. We have 
checked that these two effects do not depend on truncation effects by locally 
performing extended calculations with iVo = 14. 



4 5-type ellipsoids and superdeformed states in 50 Gd 

We now wish to present some preliminary results obtained within the gener- 
alized routhian formalism. We use the velocity field of eq. (12) which couples 
the global rotation of the nucleus with an intrinsic vortical field and is known 
after Riemann and Chandrasekhar [12] as the S-type ellipsoid approximation. 
In this particular case the expression (3) of the routhian can be developed 
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HQ. (MeV) 

Fig. 3. Single-particle routhian spectra in the 150 Gd yrast band for both protons and 
neutrons. The convention used for the (parity,signature) representation is the fol- 
lowing: (+,+) in full lines, (+,-) in dashed lines, (-,+) in dash-dotted lines and(-,-) 
in dotted lines. 



through straightforward calculations as the doubly constrained operator 

R — H — ft(Li + si) - uKi, (44) 

where K x is the first component of the Kelvin circulation operator which writes 
in cartesian coordinates 

K x = -ih \ qx 2 -^— - - XsJ—) , (45) 
\ ox 3 q dx 2 J 

q being the axis ratio as/a 2 defined in Section 2.2. This definition corresponds 
to a double stretching of the angular momentum operator in both positions 
and momenta. 

The introduction of the S-type velocity field in the context of nuclear physics 
has been already discussed by two of the authors and I. N. Mikhailov in 
Refs. [6,13]. In particular it has been stated that the Kelvin circulation de- 
fined above, which turns out to be the conjugate variable associated with the 
angle whose time derivative is called u , should be quantized as is the angular 
momentum associated to the rotation angular velocity fl This quantization 
effect has been proposed as a tentative explanation of the striking 2 h stagger- 
ing phenomenon observed in some superdeformed bands in the A = 150 mass 
region [20] and less clearly in the A = 130 mass region for Cerium isotopes 
[21]. 
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Fig. 4. Dynamical moment of inertia of / = 50 fi solutions of the nuclei 50 Gd as 
a function of the Kelvin circulation J. The dashed line represent the yrast value of 
the moment of inertia. In all the calculations reported here we have used the basis 
parameters optimized for the yrast I = 50 % solution. 



The expectation value of the Kelvin circulation is easily obtained from the 
expression (41) of the total angular momentum by performing the stretch- 
ing in coordinates and momenta and removing the spin degrees of freedom 
contribution, as 



27]- roo roo 

(K 1 /h) — — - / rdr / dz 
A Jo Jo 



-J {1) 

q Z 



The quantization rule for the Kelvin circulation writes 

(Ki/h) 2 — J(J + 1), 



(46) 



(47) 



where J must be an even integer due to the C2-symmetry character of the 
solution. 



Upon varying the Kelvin circulation from its yrast value (which corresponds 
to a vanishing uj value), it is possible to obtain a wide class of collective 
flows, including irrotational flow, shear modes and others, corresponding to 
various moments of inertia differing from the one obtained on the yrast line. 
Such effects being also observed whenever pairing correlations are introduced, 
we can thus simulate the consequences of pairing forces in terms of current 
patterns by shifting the Kelvin circulation from its yrast value through the 
second dynamical constraint —ujK\. 

To illustrate this point, we present in Fig. 4 the variation of the dynamical 
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moment of inertia as a function of J for 150 Gd solutions with / = 50 % (cor- 
responding to hQ = 0.544 MeV on the yrast line). The moment of inertia is 
seen to fit rather nicely on a parabola and can vary above and below the yrast 
value as J is varied. More precisely, classifying the solutions by the rigidity r 
of their flow patterns defined [22] as 

r = l + ^(?+l/?), (48) 

where r = corresponds to purely irrotational flow and r = 1 corresponds 
to global rotation, we can roughly associate the higher moments of inertia 
with r < 1 and the lower ones with r > 1. It is well known, as shown for 
instance by Durand et al. in Ref. [23] that the introduction of pairing forces 
leads to irrotational- like current lines. It is then very promising to notice that 
in the case of 150 Gd both pairing correlations (as shown in Ref. [19]) and 
r < 1 S-type solutions correspond to a higher moment of inertia than the 
yrast (Hartree-Fock) value. 

We will not enter much into details here, leaving for a future publication [24] 
a more complete discussion on this subject. Let us merely mention that the 
Kelvin circulation value required to obtain a moment of inertia close to the 
experimental value which is around 90 h 2 /MeV in the ~ 0.55 MeV region 
is J ~ 25 H. Such a value corresponds to a ratio J/I = 2 which is, in the 
hypothesis of Ref. [13], a condition required for a 2 h staggering to appear. 



5 Conclusions 

The aim of this paper was to present a new formalism for solving the Skyrme- 
Hartree-Fock problem in the case of a generalized routhian which breaks both 
time reversal and axial symmetries. Even though such a kind of formalism was 
already existing as we recalled in the introduction, previous works have only 
described rotating solutions (i.e. being restricted to the single usual cranking 
case). More importantly, the choice here made to use axially symmetric basis 
states, has been shown to reduce the computation time with respect to the 
usual choice of triaxial basis states. The results of our formalism in the crank- 
ing case have been tested against the results available in the literature in the 
A = 150 mass region. They are in quite satisfactory agreement for usually 
considered observables in superdeformed nuclei. 

We have also performed preliminary calculations within the S-type ellipsoid 
approximation. In Ref. [6], we had tackled the problem in the limiting case 
where the mean field is mocked-up by an harmonic potential. Here we started 
to lift up this approximation. However, this study is by no means complete 
and will be investigated in more details in forthcoming publications. It al- 
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ready gives us some hints, however, about the strong connection existing, as 
predicted, between the pairing correlations and the intrinsic vortical collective 
modes. For this purpose, it is clearly required that our formalism be extended 
to include pairing correlations even though some important results can cer- 
tainly be obtained already within the Hartree-Fock approximation. This work 
is in progress, as we are currently working in upgrading our formalism to an 
Hartree-Fock-Bogoliubov approach, using a zero-range pairing force, treat- 
ing moreover the complicated problem of projection on good particle number 
in the approximate and widely used Lipkin-Nogami scheme (see e.g. [25] for 
the implementation of this scheme in the framework of self-consistent cal- 
culations). Such a HFB formalism will also be instrumental in testing the 
hypothesis of Ref. [13] on the appearance of a 2 h staggering in some superde- 
formed bands. Indeed, since the moment of inertia renormalization induced by 
the pairing force is accompanied by a more pronounced irrotational character 
of the current patterns, the Kelvin circulation spectra along the yrast line 
should be very different in the HF and HFB cases. One should not therefore 
consider the J/ 1 value in the pure Hartree-Fock the final word and 

rather consider it as a good starting point for a more comprehensive study of 
high spin nuclear dynamics. 
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A Routhians 

A.l Skyrme energy functional and local densities 

The energy functional in the case of time-reversal symmetry breaking writes 
as [8] 

ft(r) = ftidn(r) + ttvoi(r-) + ^ so (r) + X odd (r) + tt cou i(r), (A.l) 
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where these different terms are given below, omitting their r dependence. The 
kinetic energy part writes 



taking approximately into account the center of mass energy correction. The 
rest is given in terms of the force parameters Bi (i — 1, 11) given in Ref. [26] 



Mvoi = B lP 2 + B 2 J2p 2 q + B 7P a+2 + B 8 p a J2 p\ 
i i 

+B 3 pr + 5 4 E Pfo + B bP V 2 p + 5 6 E P^Pg, (A.3) 



(A.4) 



Modd = B 9 ( V x p ■ j + V x p q ■ j q J - B 3 j 2 ~B 4 J2 f q 
V q J q 

+B 12 p a p 2 + B 13 p a J2pI + B wP 2 + B u J2 Pi (A.5) 

in which some small contributions have been neglected as was discussed in 
[26] and 

using the Slater approximation [27] for the exchange term whose accuracy has 
been checked in Ref [28]. These expressions make use of the following local 
densities 



p ff = (A.7) 

k 

^ = E^i-^, (A.8) 
k 

V-J 5 = -i^V^-Vx^ (A.9) 

k 

i,=^E $ i v$ *- $ * v$ ;» (a.io) 

P ff = E $ *^*. (A-ll) 
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where all the sums run over the occupied states for a given charge state q. 
When no charge index is present, we refer to the total density like e.g. 



(A.12) 



A. 2 Form factors of the one-body routhian 



The form factors appearing in eq. (7) are expressed as follows: 



U q = 2(B lP + B 2Pq ) + (B 3 t + B A r q ) + 2(B 5 Ap + B 6 A Pq ) 

+ (2 + a)B 7 p l+a + Bsp^ 1 (2pp q + «E^) + B ^ ' 3 + V ' J <?) 



+ap a - 1 B l2 p 2 + B13 £ pj ) + S qp Koui - e 2 



7T 



Pp 



1/3N 



(A.13) 



2 ?77 



(A.14) 



= + 2{B 3 j + B 4 j ) - 5 9 (V xp + Vxp), 



(A.15) 



= hi - B 9 (V x j + V x - 2(B 10 p + B n p 9 ) 
-2p a (B 12 p + B 13 p q ), 



(A.16) 



5, 



v; so = -^(p + p,), 



(A.17) 



and 



T/dir _ p 2 /* Pp( r ') j„/ 
^coul — e 



rdr' 



r — r 



(A.18) 



B Expression in the space coordinates of densities and energy form 
factors 
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B.l Densities calculation 



The calculation of the densities and routhian matrix elements requires the 
knowledge of the derivatives of those functions with respect to the three cylin- 
drical coordinates. The derivation with respect to 6 is straightforward from 
equation (19), whereas the other derivatives are given with obvious notations 
by 



dM{r) 



/^i c -(g 2 +») 

2tt 



1/2 



1/2 



Y,(±rcy iAe v m - 1)/2 H n A A }, (B.l) 



using as was done in Ref. [9] the definitions 



(B-3) 



L^in) = 2(n r + l)4 A J + i(^) " (2^ + |A| + 2 - V )L n M- (B.4) 



From the expressions given by Vautherin (eqs. (5.5) of Ref. [9]) for the r 
and V 2 p densities in the case of time-reversal symmetry and the calculation 
of p performed in sect. 2.2 we can easily obtain in our case the following 
expressions: 



r = M e -« 2 + .) £ p^^MA'l)/*-^, cos [(A - A')9] 

{ V fiH n ,H n ,LMLM + fi_H n .H n . [LWD£ + AA'4 A J<'] } , (B.5) 



V 2 p = 2r, 



2(3 z (3j 



7Y 



■e- (e+v) E feV' A|+|A1)/2 ^ cos[(A - A')0] 



{01 [e - 2(n z + I)] + /5i [77 - 2(2n r + |A| + 1)]} 

n. 2 ii> z n r n' 



(B.6) 



Using the expression (3.7) of Ref. [9] for the cylindrical components of the 
operator V x cr, we obtain here the divergence of the spin-orbit density as 
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v • j = 2$>{ - W^w} 

{s (e- i6 d r *fd z n - e i6 d r ^*d z *t)} , (B.7) 

where 3?{x} resp.) represents the real part (imaginary part resp.) of x. 

This expression may be worked out with the help of the equations (19), (B.l) 
and (B.2) as 

V J= M e -(5 2 +.) £ ^(IM+IA'D/2 

^A'H n H n ,Ll A }L l n f cos[(A - A')0] (l + (-)-+<) (B.8) 
-sp g p±V- 1/2 H n ,H n , ALMlM cos[(A + A' + 1)0] ((-)< - (-)»•) 
+sP z p ± r]- 1/2 H nz H n , LWlM cos[(A + A' + 1)5] ((-)< - (-)»•) . 

Now introducing the definition 

2p+1 J 1 if A - A' is an odd integer, 
I (J it A — A is an even integer, 

we finally write this density in the following way: 



V J= ^le-ie+r,) J- (|A| + |A'|-D/2 

A'^V- 1/2 H nz H n ,X A }L% l 5t AI cos[(A - A')6>] (B.10) 
+s{-Y'*(3 z H nz H n ,l)^ (ZI A J - A4 A J) cos[(A + A' + 1)61]. 

From the definition (A. 10) of the current density, one can write the r compo- 
nent of this density like 

3r = £3 {W + W} , (B.11) 

which can be readily developed as 



sin[(A' - A)9] (l - (-)"«+"'.) , (B.12) 
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since the sine is an odd function. This expression can be simplified as 



Jr = 



^e-(^) £ P ^ m+m - m ^ sm[(A - A') 



71 



H nz H n ,L^L[f. (B.13) 
In the same way, the two other components of this density are expressed as 



36 = ^e-« 2 +") J2 P^ A]+lA ' ] - 1)/2 S 2 rt cos[(A - m 

A'H nM H n ,LmM (B.14) 



: _M± e - {e+ v) £ p w V |A|+|A,| - 1)/2 e + A< sm[(A - A')0] 
^^4 A j4f- (B.15) 



To calculate the spin-vector density, we will use the cylindrical coordinates of 
the operator cr which are given in terms of the Pauli matrices by 



o r = \ [a + e~ w + (T_e if? ) , 
o 9 = -'\ (a + e- ie - a_e ie ) , 

the z component being evident. We can then write 



(B.16) 



p r = 2^{s*r*- k e- w }, (B.17) 

k 

Pe = 2j2%{s*r* k e- ie }, (B.18) 

k 

k 

It is then straightforward to obtain 



Pr = ^e-(^) £ S (-)<^V |A|+|A1)/2 ^ cos[(A + A' + 1)6] 

E nz H n ,L\^,\ (B.20) 

since the sum of the pp! and p'p terms cancel when n z + n' z is an odd integer. 
The 9 component is simply given by 
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Pe = -^le-(^) £ s{-)^ P ^ m+mi ^l P -A' H(A + A' + 1)*] 

H n ,H n ,LWLM (B.21) 



and since equations (B.19) and (20) only differ by a minus sign between the 
two terms of the right-hand side we have 



Pz = M e -(^) £ p^' A|+ ' A '' )/2 ^ cos[(A - A')9] 
ii,n' 



iiz fi-z n r n r 



(B.22) 



The components of the curl of j and p which also enter the one-body routhian 
can be obtained upon mixing partial derivatives in the following way: 



cm\ r j = d e j z - d z j e , 
cm\ e j = d z j r - d r j z , 

curl^j = d r j e - d e j r + -j e . 

r 



(B.23) 
(B.24) 

(B.25) 



B.2 Calculation of the Coulomb potential 



It is well known and used in Ref. [9] that the Coulomb potential of equation 
(A. 18) can be written after two integrations by part as 



Ko^) = -y / \r-r'\V 2 Pp (r>) dr'. 



Developing 



r — r 



f\2 



z - z'Y + (r + r ) 



1/2 



1 - k 2 cos 2 



'6-6'^ 



1/2 



with 



k 2 = 



{z - z') 2 + (r + r') 2 ' 
and writing as suggested by Table 1 

V 2 p p (r') = J2 D 2n (r',z') cos2n#', 

n>0 



(B.26) 



(B.27) 



(B.28) 



(B.29) 



the Coulomb potential can be developed in the following way 
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POO POO r -i 1 /9 

V^(r) = £ - £ / r'dr' / dz' (s - z'f + (r + rf D 2 "(r', 

Z n>0 J ° J -°° 



r 2ir 

/ d0' 

JO 



nl/2 



1 - /c COS 



cos2n9'. 



[B.30) 



In terms of u — (9 — 6')/2, the angular part of the above integral rewrites as 

I 2n = 2 J\u (l - A; 2 cos 2 w) 1/2 cos(4nw- 2nd). (B.31) 

Developing the last cosine of this expression, we can factorize the 9 dependence 
out of the integral and obtain 

I 2n = 4 cos 2n9 du (l - k 2 cos 2 u) 1/2 cos 4nu, (B.32) 

since the other term is vanishing. We can then write the Coulomb potential 
as a Fourier series in the 9 variable: 



(r) = 4e 2 E cos 2n9 / r'dr' / dz' 



n>0 



tt/2 



[z - z') 2 + (r + r') 2 



D 2n (r', z') I ' du (l - k 2 cos 2 uj 1/2 cos 4nw, 



1/2 



^B.33) 



where we have restricted the integral over z to positive values using parity 
properties. The r' and z' integrals are then performed numerically through 
10 points Gauss-Hermite and Gauss-Laguerre quadrature formulas whereas 
the u integral is performed through a 48 points Gauss-Legendre quadrature 
formula during the first iteration, the results being then stored for future use. 



C Calculations of routhian matrix elements 



Comparing our expression (32) for the scalar part of the routhian matrix 
element with Vautherin's equation (4.22) in Ref. [9], we can straightforwardly 
obtain the kinetic part of the matrix element as 



(fis\Vf ■ V|//s> = -(1 + 5°_ A ,) / dr, / d£e-« +>)r, 

Jo Jo 



-tt 2 +r?W|A| + |A'l)/2 



./ 



(|A-A'|) 



^-^^(4 a J4v' + aa'4 a J4v') 



(C.l) 



where we use, as will be done throughout this Appendix, the notation 

to represent the n-th non-vanishing component of the Fourier series of any 
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routhian form factor T (see Table 2). The calculation of the ct part of the 
matrix element can be rewritten through an integration by parts as 



{fj,s\(ct • V + V • a)|//s) = (fis\a ■ V|//s) - (fjfs\ct ■ V|//s)*, (C.2) 

which allows us to only calculate the first term on the right-hand side, the 
second term being easily deduced from it. 



(Ha • V|//s) = W ( 0i| cos[(2n + 1)5] of t+1) ^V) 

Z n>0 V r 



+ (/i| sin[(2n + 1)5] («( 2n+1) 9 r + al 2n+1) 9 2 ) |//) 



+(-)»-+■* (/x, /x') — (/i, 



(C.3) 



where the term (//, //) — * (//, //) represents a term identical to the previous 
one with p, and fx' replacing /i and //. Making use of equation (30) and of the 
following identity: 



/■27T 

Jo ' 



e ike sin m# 







17T — 



if m 7^ | A; | and m^O, 
if m = 



(C.4) 



one can notice that the \x[i! and //// contributions in equation (C.3) are iden- 
tical. The sum over n therefore reduces to 



(jjis\a ■ V|//s) = sin(|A - A'\0)(a$ A ~ A '^d r + a« A - A ' n d z ) |//) 
+(jji\cqs[\A-A'\6] a ( J J 



r 



(C.5) 



which expands as 



P OO P oo 

(fis\a ■ V|//s) =i/3jL / dr? / d£e-^V |A|+|A ' M)/2 #n.4 A| 

JO JO r 

(C.6) 



A'- A 
I A' - A| 



+e*r A ' l) AU,4f 



The spin part of the matrix element involves a spin diagonal term and two 
spin anti-diagonal terms. Using the relations (B.16) and performing the spin 
part of the scalar product, they respectively write 
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(fis\S z a z \fi's) = - £>| cos[(2n + 1)61] Sf^V) 

1 n>0 

_(_)«.+"', ^| cos[ ( 2n + i)0] ,S( 2n+1 )|/i'), (C.7) 



(/is|S r a r |/i's> = - E s(-) < (^|e~ if? cos[(2n + 1)0] ^ 2n+1 V> 

2 n>0 

+s(-) nz (/i|e ie cos[(2n + 1)0] Sl 2n+1) \fi') , (C.8) 

and 



(H^|// S > = s(-) < (/i|e- if? sin[(2n + 1)5] Sf^V) 

2 n>0 

-s(-)" z (^|e ie [sin(2n +1)5] Si 2n+1) |//) . (C.9) 

In all these three expressions again, the two terms on the right-hand side are 
the same and we finally obtain the expression 



P OO P oo 

WS.<r\v!s)= d V dee-^^l^'l)/ 2 ^^^^'! 

Jo Jo r 

S (|A-A'D + „(_)». ( S (|A+A' + 1|) 



A + A' + l (|A+A'+l|r 

lA + A' + lf* 



(CIO) 



One can easily develop the spin-orbit term of the one-body routhian of eq. (7) 

as 



U s ° = (W so x V) ■a=(d r V so )(^z ~ d z a e ) + \d 9 V so ) (d z a r - d r a z ) 

+(d z V s °)(d r a e -^a r ). (C.ll) 

As the calculation of the matrix element is rather long, we will not enter into 
details here. It has been established in Ref. [14] using the same technique as 
we used in this Appendix that, after integrating by parts to eliminate the 
(diV so ) terms, we get: 
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{lis\U so l/i's) 



- s (-)-^ 1/2 (i + W + ii)Kf +A ' +1|) K^4f' 
(4 A J - a4 a J) - ^^4 A J(4f - A'4f ')] } 




(C.12) 



References 

[I] H. Flocard, P. Quentin, A. K. Kerman, and D. Vautherin, Nucl. Phys. A 203 
(1973) 433. 

[2] J. L. Egido and L. M. Robledo, Phys. Rev. Lett. 70 (1993) 2876. 

[3] B. Girod, J.-P. Delaroche, J.-F. Berger, and J. Libert, Phys. Lett. B 325 (1994) 
1. 

[4] J. Dobaczewski and J. Dudek, Phys. Rev. C 52 (1995) 1827. 

[5] J. Bartel, P. Quentin, M. Brack, C. Guet, and H. B. Haakansson, Nucl. Phys. 
A 386 (1982) 79. 

[6] I. N. Mikhailov, P. Quentin, and D. Samsoen, Nucl. Phys. A 627 (1997) 259. 

[7] D. Samsoen, P. Quentin, and I. N. Mikhailov, in preparation. 

[8] Y. M. Engel, D. M. Brink, K. Goeke, S. J. Krieger, and D. Vautherin, Nucl. 
Phys. A 249 (1975) 215. 

[9] D. Vautherin, Phys. Rev. C 7 (1973) 296. 

[10] M. Abramowitz and I. A. Stegun. Handbook of Mathematical Functions. Dover, 
New- York, 1965. 

[II] A. L. Goodman, Nucl. Phys. A 230 (1974) 466. 

[12] S. Chandrasekhar. Ellipsoidal Figures of Equilibrium. Dover, New- York, 1987. 

[13] I. N. Mikhailov and P. Quentin, Phys. Rev. Lett. 74 (1995) 3336. 

[14] D. Samsoen. PhD thesis, Universit Bordeaux-I, 1997, unpublished. 

[15] B. Gall, P. Bonche, J. Dobaczewski, H. Flocard, and P. H. Heenen, Z. Phys. A 
348 (1994) 183. 

[16] J. Dobaczewski and J. Dudek, Comp. Phys. Comm. 102 (1997) 166 ; 183 
[17] J. Libert, private communication. 
[18] J. Meyer, private communication. 



31 



[19] P. Bonche, H. Flocard, and P. H. Heenen, Nucl. Phys. A 598 (1996) 169. 

[20] S. Flibotte et al., Phys. Rev. Lett. 71 (1993) 4229 ; 
D. S. Haslip et al., Phys. Rev. Lett. 78 (1997) 3447. 

[21] A. T. Semple et al., Phys. Rev. Lett. 76 (1996) 3671. 

[22] G. Rosensteel, Phys. Rev. C 46 (1992) 1818. 

[23] M. Durand, P. Schuck, and J. Kunz, Nucl. Phys. A 439 (1985) 263. 

[24] D. Samsoen, P. Quentin, and I. N. Mikhailov, Phys. Rev. C, in press. 

[25] P. Quentin, N. Redon, J. Meyer, and M. Meyer, Phys. Rev. C 43 (1991) 361. 

[26] P. Bonche, H. Flocard, and P. H. Heenen, Nucl. Phys. A 467 (1987) 115. 

[27] P. Gombas, Ann. Phys. (Leipzig) 10 (1952) 253. 

[28] C. Titin-Schnaider and P. Quentin, Phys. Lett. B 49 (1974) 397. 



32 



